#include <cuda.h>
#include <cuda_runtime.h>
#include <driver_functions.h>

#include "IVCalculator.h"

#ifndef Pi
#define Pi 3.141592653589793238462643
#endif

#define NUM_EXP_TIME 24

extern int procId;

//struct option;
struct GlobalConstants {

  int numOptions;
  int cudaDeviceYear;
  int cudaDeviceMonth;
  int cudaDeviceDay;
  float cudaDeviceSpot;
  float cudaDeviceRisk;
  option * cudaDeviceOptions;
};

// constant memory for optimization
__constant__ GlobalConstants cuConstIVCParams;

// device memory for results
float * cudaDeviceVolatility;


__device__ __inline__  float
CND( float X )
{

  float L, K, w ;

  float const a1 = 0.31938153, a2 = -0.356563782, a3 = 1.781477937;
  float const a4 = -1.821255978, a5 = 1.330274429;

  L = fabs(X);
  K = 1.0 / (1.0 + 0.2316419 * L);
  w = 1.0 - 1.0 / sqrt(2 * Pi) * exp(-L *L / 2) * (a1 * K + a2 * K *K + a3 * pow(K,3) + a4 * pow(K,4) + a5 * pow(K,5));

  if (X < 0 ){
    w= 1.0 - w;
  }
  return w;
}

// The Black and Scholes (1973) Stock option formula
__device__ __inline__ float
BlackScholes(char CallPutFlag, float S, float X, float T, float r, float v)
{
  float d1, d2;

  d1=(log(S/X)+(r+v*v/2)*T)/(v*sqrt(T));
  d2=d1-v*sqrt(T);

  if(CallPutFlag == 'c')
    return S *CND(d1)-X * exp(-r*T)*CND(d2);
  else
    return X * exp(-r * T) * CND(-d2) - S * CND(-d1);
}

__device__ int
julian(int year, int month, int day) {
  int a = (14 - month) / 12;
  int y = year + 4800 - a;
  int m = month + 12 * a - 3;
  if (year > 1582 || (year == 1582 && month > 10) || (year == 1582 && month == 10 && day >= 15))
    return day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
  else
    return day + (153 * m + 2) / 5 + 365 * y + y / 4 - 32083;
}

__device__ int
dateDiff(option *opt){
  int year1 = opt->year;
  int month1 = opt->month;
  int day1 = opt->day;

  int year2 = cuConstIVCParams.cudaDeviceYear;
  int month2 = cuConstIVCParams.cudaDeviceMonth;
  int day2 = cuConstIVCParams.cudaDeviceDay;

  int julian_day1, julian_day2;
  int difference;

  julian_day1 = julian(year1, month1, day1);
  julian_day2 = julian(year2, month2, day2);

  difference = abs(julian_day2 - julian_day1);

  return difference;
}
__device__ int hash(int year, int month)
{
    int index;

    int curr_year = cuConstIVCParams.cudaDeviceYear;
    int curr_month = cuConstIVCParams.cudaDeviceMonth;

    int year_diff = year - curr_year;
    int month_diff = month - curr_month;

    if (year_diff == 0) {
        index = month_diff;
    }
    else if (year_diff == 1) {
        index = 12 - curr_month + month;
    }
    else {
        index = 12 - curr_month + 12 + month;
    }

    return index;
}

__global__ void
getIV(float * resultIV)
{
  int index = blockDim.x * blockIdx.x + threadIdx.x;
  if(index >= cuConstIVCParams.numOptions)
      return;
  option opt = cuConstIVCParams.cudaDeviceOptions[index];
  for(int i = 0;i < NUM_EXP_TIME;i++)
      resultIV[i] = 0.f;
  __syncthreads();
  float s = cuConstIVCParams.cudaDeviceSpot;
  float r = cuConstIVCParams.cudaDeviceRisk;
  float t = (float)dateDiff(&opt) / 365.f;
  float k = opt.strike;
  float p = (opt.ask + opt.bid) / 2.f;

  float cpTest = 0;
  float IV = 1.f;

  float upper = 1.f;
  float lower = 0.f;
  float range = fabs(lower - upper);

  int counter = 0;

  while(1)
    {
      cpTest = BlackScholes('c', s, k, t, r, IV);

      if(cpTest > p)
        {
          upper = IV;
          IV = (lower + upper) / 2;
        }
      else
        {
          lower = IV;
          IV = (lower + upper) / 2;
        }
      range = fabs(lower - upper);
      if( range < 0.0001)
        break;

      counter++;
    }
    resultIV[hash(opt.year,opt.month)] = IV;
}
/////////////////////////////////////////////////////
extern "C" 
IVCalculator::IVCalculator()
{
  numOptions = 0;
  cudaHostSpot = 0.f;
  cudaHostRisk = 0.f;

  cudaHostOptions = NULL;
  cudaHostVolatility = NULL;
}
extern "C" 
IVCalculator::~IVCalculator()
{

  delete [] cudaHostOptions;
  delete [] cudaHostVolatility;

  cudaFree(cudaDeviceOptions);
  cudaFree(cudaDeviceVolatility);
}

extern "C" 
void
IVCalculator::loadOptionData(std::vector<option> * opts,int n, float s, float r,
			     int cy, int cm, int cd)
{

  // initialize the settings
  numOptions = n;

  cudaHostYear = cy;
  cudaHostMonth = cm;
  cudaHostDay = cd;

  cudaHostSpot = s;
  cudaHostRisk = r;

  cudaHostOptions = new option[numOptions];
  cudaHostVolatility = new float[NUM_EXP_TIME];
  for (int i = 0; i < NUM_EXP_TIME; i++)
      cudaHostVolatility[i] = 0;
  for (int i = 0; i < numOptions; i++)
      cudaHostOptions[i] = opts->at(i);
}
extern "C" 
void
IVCalculator::setup()
{
  cudaMalloc(&cudaDeviceOptions, sizeof(option) * numOptions);
  cudaMalloc(&cudaDeviceVolatility, sizeof(float) * NUM_EXP_TIME);

  cudaMemcpy(cudaDeviceOptions, cudaHostOptions, sizeof(option) * numOptions, cudaMemcpyHostToDevice);

  printf("cudaHostYear : %u | cudaHostMonth : %u | cudaHostDay : %u\n", cudaHostYear, cudaHostMonth, cudaHostDay);


  GlobalConstants params;
  params.numOptions = numOptions;
  params.cudaDeviceYear = cudaHostYear;
  params.cudaDeviceMonth = cudaHostMonth;
  params.cudaDeviceDay = cudaHostDay;
  params.cudaDeviceSpot = cudaHostSpot;
  params.cudaDeviceRisk = cudaHostRisk;
  params.cudaDeviceOptions = cudaDeviceOptions;
  cudaMemcpyToSymbol(cuConstIVCParams, &params, sizeof(GlobalConstants));
  cudaThreadSynchronize();
}
extern "C" 
float *
IVCalculator::getResultIV()
{

  cudaMemcpy(cudaHostVolatility, cudaDeviceVolatility, sizeof(float) * NUM_EXP_TIME, cudaMemcpyDeviceToHost);
  return cudaHostVolatility;
}

extern "C" void IVCalculator::calculateIV()
{

  dim3 blockDim(32, 1);

  dim3 gridDim((numOptions/32 * blockDim.x)/ blockDim.x+1);
  
  printf("numOptions : %d\n", numOptions);
  getIV<<<gridDim, blockDim>>>(cudaDeviceVolatility);
}

